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A series of numerical simulations of thermal convection of Boussinesq fluid with infinite Prandtl 
number, with Rayleigh number 10 7 , and with the strongly temperature- and depth- dependent vis- 
cosity in a three-dimensional spherical shell is carried out to study the mantle convection of single- 
plate terrestrial planets like Venus or Mars without an Earth-like plate tectonics. The strongly 
temperature-dependent viscosity (the viscosity contrast across the shell is > 10 5 ) make the con- 
vection under stagnant-lid short- wavelength structures. Numerous, cylindrical upwelling plumes are 
developed because of the secondary downwelling plumes arising from the bottom of lid. This convec- 
tion pattern is inconsistent with that inferred from the geodesic observation of the Venus or Mars. 
Additional effect of the stratified viscosity at the upper/lower mantle (the viscosity contrast is varied 
from 30 to 300) are investigated. It is found that the combination of the strongly temperature- and 
depth-dependent viscosity causes long-wavelength structures of convection in which the spherical 
harmonic degree £ is dominant at 1-4. The geoid anomaly calculated by the simulated convections 
shows a long-wavelength structure, which is compared with observations. The degree-one (I = 1) 
convection like the Martian mantle is realized in the wide range of viscosity contrast from 30 to 100 
when the viscosity is continuously increased with depth at the lower mantle. 

I. INTRODUCTION 

Single-plate terrestrial planets like Venus and Mars without an Earth-like plate tectonics is covered by a thick 
immobile lithosphere, or cold stiff lid. It is inferred from the geodesic observations (topography and gravity) of Venus 
[Rappaport et al, 1999; Konopliv et ai, 1999] and Mars [Smith et al., 1999a; 1999b] that the spatial structure of the 
thermal convection under the lid has relatively long-wavelength in which the spherical harmonic degree is dominant 
at £ = 2-3 or lower [e.g., Schubert et ai, 1990; 1997]. In particular, as for the Mars, it is generally accepted that 
the Martian crustal dichotomy was caused by a convection system dominated by £ = 1 [e.g., Sleep, 1994; Zhong and 
Zuber, 2001]. 

In numerical simulation of mantle convection in the three-dimensional (3-D) Cartesian box geometry with wide 
aspect ratios [Tackley, 1996a; Ratcliff et al, 1997, Trompert and Hansen, 1998] and in the spherical shell geometry 
[Rat cliff et al, 1996; 1997], it is shown that a highly viscous lid is formed when the temperature-dependent viscosity 
is included in their models with the stress-free boundary condition on the top surface. As the viscosity contrast goes 
up to 10 4 -10 5 , an immobile highly viscous layer (stagnant-lid) is formed. The convection under the stagnant-lid is 
characterized by numerous, small-scale cylindrical plumes surrounded by sheet-like downwellings [Ratcliff et al, 1996; 
1997; Reese et al, 1999]. These convection patterns with high-degree modes are apparently inconsistent with the 
observations. 

Here we explore the possibility that the low-degree convection under a stagnant-lid is induced by the depth- 
dependent viscosity due to the higher viscous lower mantle. The dynamical effects of a stratified viscosity profile 
on the mantle convection without lateral viscosity variations have been studied by the two-dimensional (2-D) or 3-D 
Cartesian [e.g., Hansen et ai, 1993; Tackley, 1996b] and by the spherical shell [Zhang and Yuen, 1995; Bunge et al., 
1996; Zhong et al., 2000b] models. Bunge et al. [1996] have shown that a modest increase in the mantle viscosity 
with depth has a remarkable effect on the convection pattern, resulting in a long-wavelength structure. However, 
another important factor for the mantle viscosity, i.e., the strong dependence on temperature, was absent in their 
models. The purpose of this paper is to investigate the combined effects of (i) depth-dependence and (ii) strong 
temperature-dependence on the viscosity in the resulting convection pattern. 

II. SIMULATION MODEL 

The mantle convection is numerically treated as a thermal convection in a 3-D spherical shell of a Boussinesq fluid 
with infinite Prandtl number heated from the bottom boundary. The aspect ratio of the spherical shell f^/fi is 0.55, 
which is a characteristic value of the terrestrial planets, where fg and f\ are the radii of the inner and outer spheres, 
respectively. Equations of mass, momentum, and energy conservation governing the mantle convection are scaled to 
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a non-dimensional form as follows [e.g., Schubert et al, 2001], 

V-v = 0, (1) 

- Vp + V • [q (Vv + Vv fr ) } + Ra r Te r = 0, (2) 

BT 

— +v-VT = V 2 T + H r , (3) 
at 

where v is the velocity vector, p pressure, T temperature, t time, and e r is the unit vector in the r-direction. The 
superscript tr indicates the tensor transpose. The Rayleigh number Ra scaled by the thickness of the spherical shell 
D is given by, 

Ra^RaJiX =™ &Aft) \ (4) 

\r\ J KVref 

where p is the density, g gravitational acceleration, a thermal expansivity, AT(= Tbot—T top ) the temperature difference 
between the bottom temperature Tf, t on the inner sphere and the top temperature T top on the outer sphere, k thermal 
diffusivity, and fj re f is the reference viscosity (see equation (6) below). The hats stand for dimensional quantities. 
The internal heating rate H scaled by the thickness of the spherical shell D is given by, 

H.jjAV.-a*, ( 5) 

where Q is the internal heating rate per unit mass, and c p is the specific heat at constant pressure. In this study, in 
order to focus on the effects of the temperature- and depth-dependent viscosity, all the material properties other than 
viscosity (such as thermal expansivity and thermal diffusivity) are assumed to be constant. The viscosity r\ depends 
on the temperature T and depth d as 

r]{T, d) = r)ref{d) exp [-E (T - T ref )] , (6) 

where rj re f (d) is the viscosity at the reference temperature T — T re f . The non-dimensional "activation parameter" E 
represents the degree of viscosity contrast between the top and bottom surfaces. The velocity boundary at the top 
and bottom surfaces of the spherical shell are given by impermeable and the stress-free conditions. The boundary 
conditions for T at the top and bottom surfaces are given by Tf, t — 1 and T top = 0. 

The basic equations (l)-(3) are solved by a second-order finite difference discretization. A kind of the overset 
(Chimera) grid system, Yin- Yang grid [Kageyama and Sato, 2004], is used for the computational grid (Figure 1). The 
Yin- Yang grid is composed of two component grids (Yin grid and Yang grid) that have exactly the same shape and 
size (Figure la). A component grid of the Yin- Yang grid is a low- latitude part of the usual latitude- longitude grid 
on the spherical polar coordinates. The Yin- Yang grid is suitable to solve the mantle convection problems because 
it automatically avoids the pole problems, i.e., the coordinate singularity and grid convergence that are inevitable in 
the usual latitude- longitude grid (Figure lb). Following the general overset grid method, the data on the boarder of 
each component grid are matched by mutual interpolation. All the basic quantities — v, p, T, and n — are spatially 
discretized and located in the same grid points (collocation grid method). The details of the Yin- Yang grid can be 
found in Kageyama and Sato [2004] . Sec our previous paper [ Yoshida and Kageyama, 2004] for its application to the 
mantle convection with detailed benchmark and validation tests. 

The grid points in each component grid are 102 x 54 x 158 (in r-, 6-, and (^-directions). Thus the total grid size 
for a whole spherical shell is 102 x 54 x 158 x 2 (for Yin and Yang grids). The convergences of the solutions were 
confirmed by changing the numerical resolution with 66 x 33 x 104 x 2. Time development of the convections are 
calculated until averaged quantities, such as Nusselt number and root-mean-square velocity, become stationary. 

III. RESULTS 



Calculations carried out in this paper are summarized in Table. 1. 
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A. Constant viscosity and only temperature-dependent viscosity convections 

Before we go into details of the combined effects of the temperature- and depth-dependent viscosity on the convec- 
tion, we study the phenomenology of convection pattern changes caused only by the temperature-dependent viscosity. 

Figure 2a shows a snapshot of the residual temperature of Case r7e0 in Table 1, in which the viscosity is constant, 
i.e., E = in equation (6). The Rayleigh number Ra is 10 7 , which is about one order of magnitude smaller than the 
value of terrestrial planets. (Later, we will show that the convective pattern is unchanged even when the Rayleigh 
number is increased 10 s in Case r8e6w2.) The thermal structure of Figure 2a shows a typical pattern observed in 
the 3-D spherical shell geometry. The convective flow is composed of narrow, cylindrical upwelling (hot) plumes 
surrounded by a network of long downwelling (cold) sheets. This structure is common for the purely bottom-heated 
convection. To analyze the spatial structure, the power spectrum by the spherical harmonics Y™ of temperature field 
is plotted in the right panels of Figure 2. The small scale structure (I > 10) is dominant in the middle depth, and 
the large scale structure (£ < 6) is dominant near the top and bottom surfaces that is associated with the thermal 
boundary layers. The radial profile of horizontally averaged temperature is shown in Figure 3a. The volume-averaged 
temperature is 0.26 in this case. As we expect, compared with the convections with high Rayleigh number and the 
strong internal heating [Bunge et al, 1996; Yoshida et al, 1999; Zhong et al, 2000b], the thermal structure of this 
purely bottom-heated convection are dominated by considerably long-wavelength structure. 

Figure 2b shows the results of Case r7e6r where the reference temperature T re f — 0.5. The activation parameter 
is taken to be E = ln(10 6 ) = 13.8155. The spectrum of temperature field of Case r7e6r (right panel of Figure 2b) 
shows that the power is concentrated around £ = 6-10 throughout the depth that is associated with convecting cells 
under the stagnant-lid. For this case, the volume averaged temperature is 0.72, which is larger than the constant 
viscosity convection (Figure 3b). 

In our previous paper [Yoshida and Kageyama, 2004], we did not report the cases when the viscosity of mantle 
materials has a strong temperature dependence. The regime of the flow state under the strong temperature-dependent 
viscosity in the spherical shell convection was examined by Ratcliff et al. [1996; 1997]. For the comparison with the 
previous works, the reference temperature T re f in equation (6) is fixed to the bottom temperature T(, ot in the followings. 
Therefore, the viscosity rj re f is now the viscosity at the bottom. The viscosity contrast across the spherical shell is 
defined by -y v = r](T top )/r](T bot ) = exp(E). 

Shown in Figure 4a is a regime diagram for convective flow pattern. Approximate regime boundaries are drawn. 
Our simulation results for Ratot = 10 6 -10 are shown in this diagram. The previous results by Ratcliff et al. [1997] 
(3-D Cartesian and spherical shell models) and Trompert and Hansen [1998] (3-D Cartesian model) are also included 
in the diagram. Our results basically support the previous results by Ratcliff et al. [1996; 1997]: The convecting 
pattern is classified into three regimes defined by Solomatov [1995] in the order of increasing j v ; the "mobile- lid" 
regime (Figure 4b); the "sluggish-lid" regime (Figure 4c); and the "stagnant- lid" regime (Figure 4d). 

The moderate viscosity contrast (7,, = 10 3 -10 4 ) produces the large-scale convection, or the sluggish-lid regime. 
In our previous paper [Yoshida and Kageyama, 2004], we showed that the convection at Rabot = 10 6 and j v = 10 4 
(Case r6e4) has a two cell pattern that consists of one downwelling and two cylindrical upwcllings (Figure 4b) [Ratcliff 
et al., 1995; 1996; 1997, Zhong et al, 2000b; Yoshida and Kageyama, 2004; Stemmer et al., 2004; McNamara and 
Zhong, 2005a]. In contrast, at Rabot = 10 7 and j v = 10 4 (Case r7e4), the convection pattern comes to have the 
degree-one pattern; the one cell structure that consists of a pair of cylindrical downwelling plume and cylindrical 
upwelling plume (Figure 4c). This indicates that the convecting structure in the sluggish-lid regime is sensible to the 
Rayleigh number. 

The convective flow pattern that belongs to the stagnant-lid regime emerges when j v > 10 5 . The stagnant-lid 
prevents the heat flux through the top boundary and leads to a small temperature difference in the mantle below the 
lid. The characteristic horizontal thermal structure has short wavelengths comparable to the thickness of the mantle 
(Figures 4d and 4e). This convective pattern in the stagnant-lid regime is also observed in the previous results in a 
3-D spherical shell geometry [e.g., Reese et al., 1999b]. This convective feature would be caused by the secondary 
downwelling plumes leaving from the base of stagnant-lid. At j v > 10 6 (Cases r7e6), the connected network of 
sheet-like downwelling reaches to the mid depth of convecting layer (Figure 4d). When 7,, is further increased to 
10 8 (Cases r7e8 and r7eA), the stagnant- lid become rather thick, and we clearly observe large, mushroom-shaped 
upwelling plumes (Figure 4e). 

B. Both temperature- and depth-dependent viscosity convections 

To investigate the transition in the convective pattern by adding the depth-dependent viscosity (a viscosity stratifi- 
cation), we investigate two kinds of viscosity profiles. First we examine cases in which the viscosity jumps at the phase 
transition boundary between the upper and lower mantle. Second we examine cases in which the viscosity smoothly 
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increases with depth in the lower mantle. The ratio of thickness between lower and upper mantle, cLl/cLu, is 3.39, 
comparable to that in Earth's mantle. Since the actual viscosity contrast in the depth of the terrestrial planets is not 
fully constrained, we take it as a parameter within a plausible range between 10 15 (« 30) < r) L /r] ref < 10 2 5 (« 300) 
[e.g., Davies and Richards, 1992; Karato, 2003] where 7]l is the viscosity of lower mantle. In six cases (Cases r7e6vl 
to r7e6v3, and Cases r7e6wl to r7e6w3), the initial condition is taken from the stationary state of Case r7e6r, shown 
in Figure 2b. The reference temperature T re f in equation (6) is fixed at 0.5. The Rayleigh number defined by r] re f is 
fixed at 10 7 . 

Shown in Figure 5 are the results of three cases (Cases r7e6vl, r7e6v2, and r7e6v3) in which the viscosity jumps 
at the upper/lower mantle boundary. Figure 5a shows a snapshot of the residual temperature of Case r7e6vl with 
Vh/Vref = 10 15 . Compared with the convection in which the viscosity depends only on the temperature (Figure 2b), 
we find that the convective flow pattern obviously has longer length scale. The thermal spectrum indicates a shift 
to smaller degrees, and the peak is located between £ = 2 and w 10. As r]L/T] re f is further increased (Cases r7e6v2 
shown in Figure 5b and r7e6v3 in Figure 5c), the thermal structure significantly shifts to lower modes. The power 
spectrum shows a concentration in £ < 6 with the peak of £ = 2-4 for r)L/r] re f = 10 20 (Figure 5b), and £ = 2-3 
for rjL/Vref = 10 2 ' 5 (Figure 5c) throughout the depth. As the amount of the viscosity jump i^L/Vref increases, the 
temperature drop in the bottom thermal boundary layer grows, which leads to the lower internal temperature of the 
mantle (Figure 6). 

To see this spatial scale change of the convection caused by the depth-dependent viscosity in more detail, we 
analyzed the time sequence of the Nusselt number, the root-mean-square velocity averaged over the entire mantle, 
and the peak mode at each depth for the Case r7e6v3 with f]hlf]ref = 10 2 ' 5 . At the initial stage of the simulation run, 
the convection is dominated by £ = 7-9 modes throughout the depth which reflects the initial condition (Figure 2b). 
As time goes on, the convective flow reaches to a saturated state (Figure 7a), and the low-degree component develops 
from the upper part to middle part of mantle; the peak mode shifts from £ — 9 to 3 there (Figure 7b). This indicates 
that the stagnant-lid is broken, and then, the convection cells are re-organized into the convection state with the low 
modes. 

To compare with the observation, we have calculated the geoid anomaly for Cases r7e6vl to r7e6v3. We followed 
the method of the calculation of the geoid anomaly described in Hager and Clayton [1989]. The physical parameters 
used in the calculation are set to those possibly relevant to Venus (Table 2). Figure 8 shows the distribution of 
calculated geoid anomaly where r\hl^ref is ( a ) 1 (i- e -, n0 viscosity stratification), (b) 10 15 , (c) 10 2 , and (d) 10 25 . 
The results are shown by the spherical harmonics modes up to £ — 24. Figure 8e shows the power spectrum for 
each case. The mode amplitude with the viscosity stratification peaks at £ = 2-4. When the stratified viscosity is 
absent (Figure 8a), £ = 5-10 modes are strong (see the arrow in Figure 8e). On the other hand, as rjL/i] re f increases 
(Figures 8c and 8d), the power spectrum peaks at £ — 2 and the higher degree components [t > 10) are remarkably 
decreased. This is consistent with the spectrum constructed from the observed geoid anomaly of the Venus [Konopliv 
et al, 1999] (Figure 8e). 

Next, we investigate the cases with smoothly increased viscosity with depth rather than the jump. In three cases 
(Cases r7e6wl, r7e6w2, and r7e6w3) shown in Figure 9, the viscosity contrast between the upper/lower mantle 
boundary and the bottom of mantle is Ar] L = (a) 10 15 , (b) 10 2 , and (c) 10 2 5 . The initial condition is again the 
state shown in Figure 2b. We see from both the residual temperature and the spectrum that the dominant power is 
concentrated on the smaller degrees in all the cases. At Ai]l — 10 15 -10 2 , the peak is located at £ = 1, or one-cell 
convection (Figures 9a and 9b). On the other hand, at At/l — 10 2 5 , the peak is located at £ = 2, or the two-cell 
convection (Figure 9c). The horizontally averaged temperature and viscosity profile is shown in Figure 10. Note 
that the viscosity contrast in the lid are almost identical among the three cases (see the arrow in the right panel 
of Figure 10). This suggests that the transition between degree-one and degree-two convection is sensitive to the 
magnitude of the viscosity stratification. We found that the patterns (degree-one or degree-two) are not affected by 
the increase of E up to 8 (Case r7e8w2) or to 10 (Case r7eAw2). This pattern is also unchanged when the internal 
heating is included (H = 20) (Case r7e6w2h), or the Rayleigh number is increased to 10 8 (Case r8e6w2). The patterns 
£ = 1 or 2 are mainly controlled by the viscosity contrast Atjl- 

IV. CONCLUSIONS AND DISCUSSION 

The convection with strongly temperature dependent viscosity under the stress-free boundary condition has short 
wavelength structure when the depth-dependent viscosity is ignored. This feature is inconsistent with the convection 
inferred from the geodesic observations on the single-plate planets like Venus and Mars. We have found that the 
combination of temperature- and depth-dependent viscosity produces the convection with the spherical harmonics 
degree £ = 1-4. The geoid anomaly calculated from these simulation data also generates large scale length, which 
is consistent with the observation. Schubert et al. [1990; 1997] have shown that convections with rigid boundary 
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condition on the top surface can lead to £ — 1-3 structures. In their model, however, the viscosity of fluid is spatially 
constant. Our finding is that, by considering more realistic viscosity profiles, the low-degree pattern can be reproduced 
in the convection model with stress-free boundary condition on the top surface. 

The previous convection models without the temperature-dependent viscosity [e.g., Hansen et al, 1993; Zhang 
and Yuen, 1995; Bunge et at, 1996] have already produced the large scale flow pattern by considering the viscosity 
stratification. This could be explained by the enhanced value of viscosity in the lower mantle. In our model with 
strongly temperature-dependent viscosity, the large scale convection seems to be realized by the change of convecting 
regime, from the stagnant-lid regime into the sluggish-lid regime, which is caused by the viscosity stratification. A 
major difference between their results and ours is that a highly viscous lid is naturally formed on the top owing to 
the inclusion of the temperature-dependent viscosity effect in this study. 

To date, several mechanisms have been proposed for the degree-one convection of the Martian mantle. For example, 
the endothermic phase transition just above core-mantle boundary in Martian mantle with the rigid boundary condi- 
tion [Harder and Christensen, 1996; Breuer et at, 1998; Harder, 1998], and a priori high- viscous lid [Harder, 2000] on 
the top surface boundary without any phase transitions. The small core, in other words, the thicker convecting shells 
of the mantle may lead to the degree-one convection in the ancient Mars [Schubert et al, 1990] and the Moon [Zhong 
et al, 2000a]. McNamara and Zhong [2005a] have recently found that the internal heating plays a role in increasing 
flow wavelength and forming the degree-one convection in convections in which the viscosity moderately depends on 
temperature. One of our findings in this paper is that the degree-one convection can be relatively easily reproduced 
when both effects of the temperature- and depth-dependence on the viscosity are taken into account. Although the 
degree-one convection appears even when the depth-dependence is absent (Figure 4c), the parameter range for this 
pattern is rather narrow; it is sensitive to the Rayleigh number. On the other hand, when the viscosity in the lower 
mantle is continuously increased with depth, the degree-one (I = 1) convection like the Martian mantle is realized in 
the wide range of viscosity contrast from 30 to 100. 

It is an interesting possibility that the transition of the convecting patterns between low-degree convective mode 
and one-degree mode took place in the planets. We have not directly observed the transition of convecting mode in 
our simulations. The physical parameters (Ra and/or E in this study) to characterize the convective pattern are fixed 
in our simulations. However, we would like to point out again a drastic difference of the convection patterns between 
relatively close conditions: the convection of degree-two at Ra — 10 6 (Figure 4b) [Yoshida and Kageyama, 2004], and 
the convection of degree-one at Ra = 10 7 (Figure 4c). This sensitive change has not been reported so far. 

Our simulation results will not be directly applied to the Earth's mantle, because the effects of the plate tectonics 
would be comparable to the effects of the depth-dependent viscosity, as proposed by Bunge and Richards [1996] 
and Bunge et al. [1996; 1998] from their model without temperature-dependent viscosity. Existence of a stationary 
continental lithosphere [Yoshida et al, 1999], a drifting continental lithosphere [Phillips and Bunge, 2005], and plate 
motion on the top surface boundary [Zhong et al, 2000b] also transform the small scale convection patterns in high 
Rayleigh number convection into the large scale convection patterns. 
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TABLE I: List of runs employed in this study 



Case Name Ra 


T re f 


E 


t)l/t)u 




i.e.* 1 


Corresponding figures 


r6e0 


10 6 


— 


In 10° 


— 


— 


— 


— 




r6el 


10 6 


1.0 


lnlO 1 


— 


— 


r6e0 


— 




r6e2 


10 6 


1.0 


lnlO 2 


— 


— 


r6el 


— 




r6e3 


10 6 


1.0 


lnlO 3 


— 


— 


r6e2 


— 




r6e4 


10 6 


1.0 


lnlO 4 


— 


— 


r6e3 


Fig. 


4b 


r7e0 


10 7 


— 


In 10° 


— 


— 


r6e0 


Figs. 


2a and 3 


r7el 


10 7 


1.0 


lnlO 1 


— 


— 


r7e0 


— 




r7e2 


10 7 


1.0 


lnlO 2 


— 


— 


r7el 


- 




r7e3 


10 7 


1.0 


lnlO 3 


— 


- 


r7e2 


- 




r7e4 


10 7 


1.0 


lnlO 4 


— 


— 


r7e3 


Fig. 


4c 


r7e5 


10 7 


1.0 


lnlO 5 


— 


— 


r7e4 


— 




r7e6 


10 7 


1.0 


lnlO 6 


- 


— 


r7e5 


Fig. 


4d 


r7e8 


10 7 


1.0 


lnlO 8 


— 


— 


r7e6 


Fig. 


4e 


r7eA 


10 7 


1.0 


In 10 10 


— 


— 


r7e6 


— 




r7e6r 


10 7 


0.5 


lnlO 6 


— 


— 


r7e6 


Figs. 


2b and 3 


r7e6vl 


10 7 


0.5 


lnlO 6 


lnlO 1 - 5 


— 


r7e6r 


Figs. 


5a and 6 


r7e6v2 


10 7 


0.5 


lnlO 6 


In 10 20 


— 


r7e6r 


Figs. 


5b and 6 


r7e6v3 


10 7 


0.5 


lnlO 6 


lnlO 2 ' 5 


— 


r7e6r 


Figs. 


5c and 6 


r7e6wl 


10 7 


0.5 


lnlO 6 


- 


Inl0 15 


r7e6r 


Figs. 


9a and 10 


r7e6w2 


10 7 


0.5 


lnlO 6 




lnlO 2 - 


r7e6r 


Figs. 


9b and 10 


r7e6w3 


10 7 


0.5 


lnlO 6 




lnlO 2 - 5 


r7e6r 


Figs. 


9c and 10 


r7e8w2 


10 7 


0.5 


lnlO 8 




lnlO 2 - 


r7e6w2 






r7eAw2 


10 7 


0.5 


In 10 10 




lnlO 2 - 


r7e6w2 






r8e6w2 


10 8 


0.5 


lnlO 6 




lnlO 2 - 


r7e6w2 






r7e6w2h* 2 


10 7 


0.5 


lnlO 6 




lnlO 2 - 


r7e6w2 







(*1) "I.C." indicates the Initial conditions. (*2) "Case r7e6w2h" is a case with internal heating (see text). 



TABLE II: List of parameters used in the calculation of the geoid anomaly 

Symbols Values 
outer radius n 6.052 x 10 6 m 

inner radius ro 0.55ri m 

thickness of the mantle D 0.45ri m 

density p 3.3 x 10 3 kg m~ 3 

density contrast at the top surface Apbot 2.3 x 10 3 kg m -3 

density contrast at the bottom surface Ap top 4.3 x 10 3 kg m -3 
gravity acceleration g 8.9 m s~ 2 

thermal expansivity a 1.0 x 10~ 5 K _1 

temperature difference across the mantle AT 2.0 x 10 3 K 

specific heat at constant pressure c p 1.2 x 10 3 J kg -1 K _1 

thermal diffusivity k = k/pc p 8.1 x 10~ 7 m 2 s _1 

thermal conductivity k 3.2 W m _1 K _1 

reference viscosity rj 1.5 x 10 21 Pa s 

gas constant R 8.3145 J mol -1 K _1 

gravitational constant G 6.6726 x 10~ n N m 2 kg" 2 

The values are referred with Schubert et al. [1990], Solomatov and Moresi [1996], and Turcotte and Schubert [2002]. 
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Fig.1 (Yoshida & Kageyama) 

FIG. 1: The Yin- Yang grid. Two component grids of the Yin- Yang grid are identical (the same shape and size): (a) The low 
latitude part (7r/4 < 9 < 37r/4, — 37r/4 < <p < Sn/A) of the latitude-longitude grid, (b) They partially overlap each other on 
their boarders to cover a spherical surface in pair. As it is apparent, the Yin- Yang grid has neither a coordinate singularity, 
nor grid convergence; the grid spacings are quasi-uniforms on the sphere. 
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Fig.2 (Yoshida & Kageyama) 

FIG. 2: The iso-surface of residual temperature 8T (i.e., the deviation from horizontally averaged temperature at each depth), 
and the power spectrum of the spherical harmonics of temperature field at each depth, (a) Case r7e0 with the constant viscosity 
(i.e., E = 0) convection and (b) Case r7e6r with the strongly temperature-dependent viscosity (E = lnlO 6 ) are shown. Blue 
iso-surfaces indicate ST = (a) —0.10 and (b) —0.15. Yellow indicate 8T — (a) +0.10 and (b) +0.15. The logarithmic power 
spectrum are normalized by the maximum values at each depth. White regions in maps indicate the values with lower than 
10~ 2 (see color bars). 



11 




0.0 0.2 D.4 o.e O.fl 1.0 
□veragad temparalure 



Fig,3 (Yoshida & Kageyama) 



FIG. 3: Radial profiles of horizontally averaged temperature at each depth. Two cases (a) and (b) correspond to each case 
shown in Figure 2 (Cases r7e0 and r7e6r, respectively). 




Fig, 4 (Yasfiida 6 Kageyama) 

FIG. 4: (a) The three convection regimes with varying Rayleigh number (Ra 00 t) and the viscosity contrast across the shell 
(7,,); the mobile-lid (circles), the sluggish-lid regime (triangles), and the stagnant- lid regimes (squares). Solid marks show 
our calculations. Open marks show the results from 3-D Cartesian box and spherical shell models by Ratcliff et al. [1997]. 
Gray marks show the results from 3-D Cartesian box models by Trompart and Hansen [1998]. The regime boundary (dashed 
curve) between convection regime and no-convection regime is referred with the reviews by Schubert et al. [2001]. Dashed line 
shows the approximate boundaries that separate the three convection regimes, (b)-(d) The iso-surface renderings of residual 
temperature shown in (a); (b) Ra^ot = 10 6 and 7^ = 10 4 (Case r6e4), (c) Rai, ot — 10 r and 7^ = 10 4 (Case r7e4), and (d) 
Rabat = W 7 and j v = 10 6 (Case r7e6). Blue iso-surfaces indicate (b) ST = —0.20, (c) —0.25, and (d) —0.10. Yellow indicate 
(b) ST — +0.40, (c) +0.25, and (d) +0.10. The red spheres show the bottom boundary of the mantle, (e) The temperature 
distribution on a cross section for a case where Ra oot = 10 7 and j v = 10 s (Case r7e8). 
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Fig. 5 (Yoshida & Kageyama) 

FIG. 5: The iso-surface of residual temperature (ST) and the power spectrum of the spherical harmonics of temperature field 
at each depth for the cases where r]h/r) r ef = (a) 10 1 ' 5 , (b) 10 2 , and (c) 10 2 ' 5 (Cases r7e6vl, r7e6v2, and r7e6v3, respectively). 
Blue iso-surfaces indicate (b) 5T = -0.20, (c) -0.25, and (d) -0.25. Yellow indicate (b) ST = +0.20, (c) +0.25, and (d) +0.25. 
The logarithmic power spectrum are normalized by the maximum values at each depth. White regions in maps indicate the 
values with lower than 10 -2 (see color bars). 
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Fig. 6 (Yoshida & Kageyama) 

FIG. 6: Radial profiles of the horizontally averaged temperature (left), and the horizontally averaged viscosity (right) at each 
depth. Three cases (a)-(d) correspond to each case where r)L/iJref = ( a ) 1 (i- e -, no viscosity stratification), (b) 10 1 ' 5 , (c) 10 2 , 
and (d) 10 2 ' 5 (Cases r7e6vl, r7e6v2, and r7e6v3, respectively). 
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Fig J (Yoshida & Kageyama) 

FIG. 7: The time sequence of (a) the Nusselt number (dashed line) and the root-mean-square velocity averaged over the entire 
mantle (solid line), and (b) the maximum power spectrum at each depth. The range of the spherical harmonic degrees (I) is 
analyzed up to I = 10. 
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Fig. 8 {Yoshida & Kageyama) 

FIG. 8: The contour plots of the distribution of geoid anomaly for each case where r/L/r/ref ~ (a) 1 (i.e., no viscosity 
stratification), (b) 10 1 ' 5 , (c) 10 2 , and (d) 10 2 ' 5 . The results are shown by the spherical harmonic expansion up to £ — 24. The 
spectrum are normalized by the maximum values at each degree, (e) The power spectrum of the calculated geoid anomaly for 
each case (thin colored lines) and the observed geoid anomaly from the data by Konopliv et al. [1999] (thick black line). The 
spectrum are normalized by the maximum values of all degrees. 
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Fig. 9 (Yoshida & Kageyama) 

FIG. 9: The iso-surface of residual temperature (ST) and the power spectrum of the spherical harmonics of temperature field 
at each depth for the cases where Ar/L = (a) 10 1 ' 5 , (b) 10 2 , and (c) 10 2 ' 5 (Cases r7e6wl, r7e6w2, and r7e6w3, respectively). 
Blue iso-surfaces indicate (a) ST = -0.30, (b) -0.30, and (c) -0.20. Yellow indicate (a) ST = +0.30, (b) +0.30, and (c) +0.20. 
The logarithmic power spectrum are normalized by the maximum value at each depth. White regions in maps indicate the 
values with lower than 10 -2 (see color bars). 
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Fig. 10 (Yoshida & Kage; 



FIG. 10: Radial profiles of horizontally averaged temperature (left), and the horizontally averaged viscosity (right) at each 
depth for each case where Ar/L = (a) 1 (i.e., no viscosity stratification), (b) 10 ls , (c) 10 2 , and (d) 10 2 ' 5 (Cases r7e6wl, 
r7e6w2, and r7e6w3, respectively). 



